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Abstract 



We consider the integrate-and-fire model of the cardiac pacemaker with delayed pulsatile coupling. Sufficient conditions of syn- 
chronization are obtained for identical and non-identical oscillators. 
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1. Introduction 

In paper Jjj C. Peskin proposed a model of cardiac pace- 
maker, where signal of fire arises not from an outside stimuli, 
but in the population of cells itself. Well known conjectures of 
self-synchronization were formulated. Solution of these con- 
jectures for identical oscillators tfl} 0] stimulated intensive in- 
vestigations (3]- Hi]. 

Delays arise naturally in many biological models [12]. In 
particular, they were considered in firefly models 11 1 311 as de- 
lay between stimulus and response, and in continuously cou- 
pled neuronal oscillators 01411 . Authors of [4] considered the 
phenomenon for the Mirollo and Strogatz analysis, [d]. Iden- 
tical oscillators were investigated. Two oscillators dynamics is 
discussed mathematically, and the multi-oscillatory system by 
computer simulations. It was found that the excitatory model 
of two units "can get only out-of-phase synchronization since 
in-phase synchronization proved to be not stable." In paper f^ll 
a model without leakage was discussed, that is, oscillators in- 
crease at a constant rate between moments of firing. It was 
found that a periodic solution is reached after a finite time. Con- 
sequently, research of integrate-and-fire models, which admit 
delays and fire in unison is still on the agenda. 

In paper lll5ll we have introduced a new method of investiga- 
tion of biological oscillators. The method seems to be universal 
to analyze quite identical integrate-and-fire oscillators. In the 
present article we extend it to the Peskin's model with delayed 
interaction. Conditions are found, which guarantee synchro- 
nization of the model. Our system is different than that in 04J], 
since we suppose that the pulse-coupling is instantaneous, if 
oscillators are close to each and are near threshold. In next our 
papers, we plan to consider other models, varying types of the 
delay involvement, as well as inhibitory models such that ana- 
logues of results in A and IB can be obtained. The method of 
analysis of non-identical oscillators is based on results of the- 
ory of differential equations with discontinuities at non-fixed 
moments iCll-Hl. 



2. The couple of identical oscillators 

Let us start to analyze two identical oscillators, which satisfy, 
if they do not fire, the following differential equations 



yxu 



(2.1) 



where < Jt* < l,i = 1,2. It is assumed that S,y are positive 

numbers and k = — > 1. 

r 

When x/t) — 1, then the oscillator fires, Xj(t+) = 0. The 
firing changes value of the another oscillator, x, , such that 



Xi(t+) = 0, if Xi(f) > 1 - e, 



and 



Xi(t + t+) = Xi(t + r) + e, if x,(f) < 1 - e. 
We have that 

Xi (s) = Xi(*)e" r(J "° + J e-v (s - u) Sdu 

near t. 

In what follows, assume that 



(2.2) 



(23) 



1 



k - 1 + e 
Then, from 



< e 



(2.4) 



lto(jr)|| < \\x i (t)\\e- yU - t) + J e- y{s - u) Sdu 
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(1 -e)e- yT +k(1 -e~ yT ), 

and Xj(t) < 1 - e, we obtain that Xj(s) < 1 , for all s e [t , t + t] . In 
other words oscillator x, does not achieve the threshold within 
interval [f, t + f\, if the distance of jt,(f) to threshold is more than 
e. This is important for the construction of the prototype map, 
and makes a sense of condition ( 12.31 ). 

One must emphasis that couplings of units are not only 
delayed in our model. By ( I2.21 i oscillators interact instanta- 
neously, if they are near threshold. This assumption is natu- 
ral as firing provokes another oscillator, which being close to 
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threshold "is ready" to react instantaneously. Otherwise, the 
interaction is retarded. 

Next, we shall construct the prototype map. Fix a moment 
t = when xi fires, and suppose that oscillators are not syn- 
chronized. In interval [£, £ + t] oscillator X2 moves by law 
x 2 (t) = X2(Oo' 7( '-° + Ce~ y( '~ u) Sdu, and 



xi(£ +t)= 1*2(0 - K]e yT + k. 
Denote t — tj, the firing moment of x 2 , then 

x 2 (ri) = [x 2 (C + t) + e]e-^- ( - T) + k[1 - sT*"-*^]. 
The equation X2(rf) - 1 implies that 

g-yfo-f) _ 1 ~ * 

x 2 (0 -K + e\ 

where e\ = ee yT . Since x\(rf) - k[1 - e _r ^ we have that 
1 - (Jfe(f) + d) 



(2.5) 



(2.6) 



f-fe(0 + e i) 
Introduce the following map 

1 - (y + e,) 



Ld(v, e) = k 



k-{v + €iY 



(2.7) 



(2.8) 



such that x\(rf) = Lo(x2(£)),e). If t = g is the next to 77 fir- 
ing moment of X2, then one can similarly find that X 2 (g) = 
Ld{x\{t]), e). One can see that the map Lq can be useful for 
our investigation, since it evaluates alternatively the sequence 
of values x\ and X2 at firing moments. 
Take t > so small that 



e " > e. 



(2.9) 



From (12.9b it implies that e\ < 1 . 
One can evaluate that 

L D (l-ei,e)=0. 
and the derivatives of the map in (0, 1 — ei) satisfy 



1 — K 

L' D (v, e) = K- ; — < 0, 



and 



l—ic 

L"(v, e) = 2k ; — < 



We can easily find that there is a fixed point of the map, 



V ={K-—)- y\K J --K+—, 



and 



L' D (v*,e) < -1. 
That is, fixed point v* is a repellor. 



(2.10) 



(2.11) 



(2.12) 



(2.13) 



Now, we will define an extension of Lp on [0, 1] in the fol- 
lowing way. Let 



1-e, 

0) = K . 

K-e\ 

One can see that 1— £<o><l,if 



K 



< 



1 +e 



(2.14) 



(2.15) 



In what follows, we assume that e is sufficiently small such that 
( 12.4I > implies ( 12.151 ). We set L D (0, e) = cj, and define L D (v, e) = 
0, ifl-ei <v< 1. Since Ld '■ [0, 1] — > [0, 1] is a monotonic 
continuous function and [0, 1] is an invariant set of this map, 
it is convenient for analysis by using iterations. The graph of 
this map is seen in Figure Q] To emphasize a significance of 




Figure 1: The graph of map Ld in red, fixed point v*, and stabilized trajectory 
are seen. 

this map for the present analysis, let us see how iterations of 
it can help to observe the synchronization. Fix ?o > 0, a firing 
moment, :q(fo) = 1, *i(fo+) = 0. While the couple Xi,x 2 , does 
not synchronize, there exists a sequence of moments to < t\ < 
. . . such that x\ fires at f, with even i and X2 with odd indices. 
Denote m, = x\{ti), if i is odd, and m, = x 2 (ti), if i is even. One 
can easily see that u- l+ \ = Lo{ui, e), i > 0. The pair synchronizes 
if and only if there exists j > 1 such that x\{t) + X2(t), if t < 
tj, and x\(t) = x 2 (f), for t > tj. In particular, both oscillators 
have to fire at tj. That is, inequalities 1 — e < Uj-\ < 1 are 
valid. It is possible if < Uj- 2 < L D l (l - e). In particular, 
we have that Lo(0) = oj satisfies this condition. In the same 
time, if 1 — e\ < Uj-3 < 1, then Uj- 2 = = Lo(uj-3) and 
1 - e < Uj-\ — cd < 1 again. That is, we have found that if there 
exists an integer k > such that 1 - e < L k D (v) < 1 , then the 
motion (x\(t), x 2 (f)) with x\(to+) = v, X2(to+) = 0, synchronizes 
at the k-th firing moment. Conversely, if a motion (xi(i), x 2 (t)) 



2 



synchronizes, then one can find a firing moment, to, such that 
*i(?o+) = 0,^2(^0+) = v, v e [0, 1], and a number such that 

1-6<£*(V)<1. 

Thus, the last discussion confirms that the analysis of syn- 
chronization is consistent fully with the dynamics of the intro- 
duced map Lo(v, e) on [0,1], and the map Lq can be applied 
as the main instrument of the paper. That is why, we use this 
function as a prototype map in our investigations. 

Now, by applying properties of Lb, and analyzing self- 
compositions of the map, one can easily obtain that for all k > 
functions L k D have only one fixed point, v*, and |[L^(v*, e)]'| > 
1. We skip the discussion as it is respectively simple, and re- 
quest a large place. Since all the maps L k D have one and the 
same fixed point, v*, there is not a fc-periodic motion, k > 1, of 
the map. Consequently, for arbitrary point v + v* one has a sta- 
bilized trajectory present in Figure [1] The couple synchronizes 
when L D (y, e)> 1 — e\. 

Next, we investigate the rate of synchronization. Set ao = 
L D l (u) = Q,a k+l = L D l {a k ),k = 0, 1,2, . . . (See Figure 0. 



less than T , we can conclude that the following theorem is cor- 
rect. 




Figure 2: Boundaries, a,-, of rate intervals are seen. 

Denote by Sk the region of [0,1], which points v are synchro- 
nized after exactly k— iterations of the map Lq. One can see that 
So = [1 — £1, 1], 5 1 = [ao,ci2] and Sk = (ai-i,a*+i]> if k > 3, 
is an odd positive integer, and Sk = [ak+i,ak-i), if k > 2, is an 
even positive integer. One can see that ak — > v* as k — > 00. We 
shall call S k, k > 0, the rate intervals. 

From the discussion has been made above it follows that 
there is no a finite time such that all points of the unite square 
synchronize. The closer v is to the equilibrium v* the later is the 
moment of synchronization. 

Set T = ^ In and denote by f the time needed for solu- 
tion u(t, 0, v*) of the equation u' = S — yu, to achieve threshold. 
Since all oscillators fire within an interval of length T and the 
distance between two firing moments of an oscillator are not 



Theorem 2.1. Assume that \2.4\ and \2.9\ are valid. If to > is 
a firing moment, x\(to) — l,xi(?o+) = 0, and X2(?o+) 6 S m for 
some natural number m, then the couple Xi,Xj of continuously 
coupled identical biological oscillators synchronizes within the 
time interval [to + §7\ to + Tm]. 

3. Non-identical oscillators: the general case. 

To make our investigation closer to the real world problems 
one has to consider an ensemble of non-identical oscillators. 
We will discuss the following system of equations 



x't = (S + pi) - (y + £i)Xi, 



(3.16) 



where < Xi < 1 + i = 1 , 2, . . . , n. The constants S and y are 



the same as in the last section such that k = - > 1 . Moreover, 

y 

constants p, and & are sufficiently small for k, = > 

7 0. The 



1. 



When Xj(t) = 1 + then the oscillator fires, Xj(t+) 
firing changes values of other oscillators x,-, i + j, such that 



Xi(t+) = 0, if Xi(t) > 1 - e 
and, if x,(f) < 1 — e, then 

Xi(t + t+) = Xi(t + r) + e + e,. 



(3.17) 



(3.18) 



In what follows, we call real numbers e,pi, e,-, parame- 
ters, assuming the first one is positive. Moreover, constants 
Pi,^i,^i,€i will be called parameters of perturbation. Assume 
that they are zeros to obtain the model of identical oscillators. 
We have that an exhibitory model is under discussion, that is 
e + e, > for all i. Coupling is all-to-all such that each firing 
elicits jumps in all non-firing oscillators. If several oscillators 
fire simultaneously, then other oscillators react as it just one os- 
cillator fires. In other words, any firing acts only as a signal 
which abruptly provokes a state change, the intensity of the sig- 
nal is not important, and pulse strengths are not additive. We 
have that 

Xi (s) = *i(*)e" (r+fiX '" f) + J e -<-y+Ms-»Xs +pi\3u\9) 

near t. 

If one assume that condition ( 12.4b is valid, and constants /i, 
and £ are sufficiently small such that 



k, 1- 1 + e 



< e 



-(y+£)r 



(3.20) 



then Xi(s) < 1 for all s 6 [t, t + t], if x(t) < 1 - e. 

In this section we begin with analysis of a couple of oscil- 
lators of the ensemble of n oscillators, and find that the couple 
synchronizes if parameters close to zero. Then synchronization 
of the ensemble will be proved. 

Consider the model of n non-identical oscillators given by 
relations (12. U and (12. 31 . Fix two of them, let say, x\, x r . 
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Lemma 3.1. Assume that conditions \2.4\ and \2.9\ are valid, 
to > is a firing moment such that x;(?o) = 1 + £/, x/(fo+) = 0. 
If parameters are sufficiently close to zero, and absolute values 
of parameters of perturbation are sufficiently small with respect 
to e, then the couple xi, x r synchronizes within the time interval 
[to, fo + T] if x r (to+) t [ao, «i) and within the time interval [to + 
^-f ,f + (m + 1)71, ifx r (to+) e S m ,m > 1. 

Proof. If 1 + - s - s r < x r (to) < 1 + £ r , then two oscillators 
fire simultaneously, and we have only to prove the persistence 
of the synchrony, that will be discussed later. So, fix another 
oscillator x r (t) such that < x r (to) < \ + t; r — s — s r . While 
the couple is not synchronized, there is a sequence of firing 
moments, such that < to < t% < . . . , and oscillator xi fires 
at ti, with i even, and x r fires at f, with odd For the sake of 
brevity let u t = xi(tf),i = 2j + 1,k, = x r (ti),i = 2j,j > 0. 
In what follows we shall show how to evaluate u l+ \ through 
L(uf). Consider i is even. There are k < n — 2 distinct firing 
moments of the motion x(t) in the interval (?,-, f,- + i). Denote by 
t, < 9\ < 62 < ... < 6k < fj+i, the moments of firing, when at 
least one of the coordinates of x(t) fires. We have that 

x r (6i +t) = (x r (ti + t) + e + e^Q-WrW^r-td + 
x r (0 2 +t) = (x r (6i + t) + e + £r )e- (y+ir)(ei - e,) + 



-(r+f r )(fl2-fli) 



x r (6j +t) = (x r (6j-i + t) + e + e^y^J-e^) + 



x r (t i+l ) = (x r (6k + t) + e + er ) e -(r+«fe + i-%-T) + 

K r (l - e-trHM+i-Bk-Vy (3 21) 



The moment f, + i satisfies 

1 + £. - e - e r < X r (t i+ i) < 1 + <f r , 

and continuously depends on parameters and x r (f,). 
We have also that 



(3.22) 



x l (e l +T) = Ki(l-e-W> )( - ei+T -" ) ), 
xi(02 + t) = (xi(6i + t) + e + e/)e 



-r(r+6)(82-fli) 



k,(1 -e 



(3.23) 



jc/(0/ + r) = (x,(0/-i + t) + e + e; ) e -fr + ^-%-i) . 



x,(t i+l ) = (x,(0* + t) + e + e/ )e-^')fe'-*-) + 
^(l_ e -(r+fflfe +I -ftr-r)). 

The last two formulas describe the dependence of m, + j on m,-. 
One can easily find that similar dependence can be found if i is 
odd. 



Set StifXt, £) = ki-k. One can see that <5,(0, 0) = 0. Use (EL2D 
and (13.23b to obtain 

xAf M ) = (xAft + t) + e )e-^'-'') e -^ fel -' i) + 



K {\-Q<y^r)(t M -td )+erQ - 



(r+fr)('i + i-4 



./=! 



and 



x/(fi +1 ) = (^ + 50(l-e- (T+f ' x ' M - fi) ) + 



-(7+(l)(tM-0j-T) 



(3.25) 



7=1 



Now, recall map Lq defined in the last section. We have 

4>(t M ) = (x r (ti + t) + e)z-^-<-- T) + 
k(1 - 

where f,+ i satisfies 

<f>(? M ) = 1, 

and 

^F i+1 )=41-e- 7(6+1 - , '°). 



(3.26) 



(3.27) 



(3.28) 



By applying the definition of Lo one can see that Lo{u,) = 
<A(f i+ i)- 

Assume, without loss of generality, that F !+ i < f, + i. Then one 
has that 



<f>(t M ) - x r (t i+l ) = 1 - x r (t i+ i) = 
0>i(e,£ r ,f r , S r ,r), 

where 

0> 1 (e,e r ,^,«y n T) = Wl-e-^ 1 ^)- 
(jc r (?i + t) + e ) e -rffi + i-«i)]( e -^ffi+i-« - 1) - e ,. e 



(3.29) 



-7(fi+l-»i) e -fr(5+l-«<). 



(e + e r ) £ e -(rK.)fe-e-T) _ tfr(1 _ e -(r + f,)(fc-« )) 

and the last expression tends to zero as all of its arguments tend 
to zero. Next, by utilizing ( 13.22b and ( |3.29b we have that - 
f,+ i < ( &2(e, £ r , (r, S r ), where 



2 (e, e r , <J r , t) = 



\^ r \ + e + \e r \ + ^i(e,e r ,^,S r ,T) 



S -\fx r \-y-\Zr\ 

Now, by applying the last inequality, ( 13.25b and ( 13.281 ) one can 
see that 

\L D (ui) - Ki(ui)\ = \x t (t i+ i) - ijj{t i+ i)\ < \xi(t i+ i) - xi(t i+ i)\+ 
\x0 i+l ) - iftt M )\ < «D 2 (5 + \m\ +7+ \m + *i- 

That is, difference Lo(uj, e) - Ui+\ can be made arbitrarily small 
if the parameters are sufficiently close to zero. Moreover, we 
should assume smallness of absolute values of the parameters 
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of perturbation with respect to e, to satisfy ( 13.221 ). This con- 
vergence is uniform with respect to uq. We can also vary the 
number of points 6, and their location in the intervals (tj, tj+i) 
between and n - 1. The convergence is indifferent with respect 
to these variations, too. 

Consider L' D (uo, e). It is true that L"^{uo, e) e [1 - e, 1]. As- 
sume, without lost of generality, that m is an even number. 
Since Lq is a continuous function, we can find recurrently, by 
applying the following sequence of inequalities \uj-L' D (uo, e)\ < 
\ui-L D (ui- U e)\ + \L D (ui-i,e)-L D (L' D l (uo, e))\,i = 1,2, . . ., that 
either 1 +% r -e-e r < u m < 1 +% r or 1 e— e; < u m+ \ < 1 
if the parameters are sufficiently small. From the notation it im- 
plies that each of the last two inequalities bring the couple to 
synchronization. 

Since each of the iterations of Lo is done within interval with 
length not more than T, we obtain now that the couple xi, x r is 
synchronized not later than f = to + (m + l)T. 

We have found that oscillators x\ and x r fire in unison at some 
moment t — 6. Next, we show that they will save the state, being 
different. To find conditions for this, let us denote by t > the 
next moment of firing of the couple. Let say, x r fires at this 
moment. Thus, we have that xi(8+) - x r {6+) = 0. Then X[(t) = 
x r (t), 8 < t < t. It is clear that to satisfy Xi(t+) = x, (r+) = 0, we 
need 1 + £ r - e - e r < x;(r). By applying formula (13.221 1 again, 
this time with ?,■ = 0, ? J+ i = r, one can easily obtain that the 
inequality is correct if parameters are close to zero and absolute 
values of the parameters of perturbation are small with respect 
to e. Thus, one can conclude that if a couple of oscillators is 
synchronized at some moment of time than it continues to fire 
in unison for ever. The lemma is proved. 

Let us extend the result of the last Lemma for the whole en- 
semble. 



Theorem 3.1. Assume that M.4\ and \2.9\ are valid , to > is 
a firing moment such that Xj(to) = 1 + £j, xj(to+) — 0. If the 
parameters are sufficiently close to zero, and absolute values of 
parameters of perturbation are sufficiently small with respect 
to e, then the motion x(i) of the system synchronizes within the 
time interval [to, to + T], if x,(fo+) & [ao,a\),i + j, and within 
the time interval [to + max, ^ t '' 1 f j () + (max^fcj + 1)T], if there 
exist x s (to+) 6 [ao, a\)for some s + j and x,(fo+) £ S k,, i + f 

Proof. Consider the collection of couples (x, , xf), i + j. Each 
of these pairs synchronizes by the last Lemma within interval 
[fo + max ' * 2 J - 1 1, to + (maxtfj ki + 1)7"]. The theorem is proved. 

Let us introduce a more general system of oscillators such 
that Theorem [XT] is still true. 



and, if x,(f) < 1 — e, then 

k 

Xi(t + t+) = xt(t + r) + e + ^ 



(3.31) 



Consider a system of n oscillators given such that if i-th os- 
cillator does not fire or jump up, it satisfies i— th equation of 
system ( 12. 11 1. If several oscillators x- u , s = l,2,,..,k, fire such 
that Xi t (t) = 1 + <p(t, x(f), x(t - Ti s ), where \<p(t, x(t), x(t - r,)| < 
iji, i = 1,2, ...,n, and x 1; (f+) = 0, then all other oscillators 
Xi , p = k + 1 , k + 1 , . . . , n, change their coordinates by law 



Xi(t+) = 0, if Xi(t) > 1 - e 



(3.30) 



One can easily see that the last theorem is correct for the 
model just have been described, if e+ 2* =1 £i p i s > 0, for all pos- 
sible k, i p and i s , and we assume that ey are also parameters of 
perturbation. Moreover, one can easily see that initial functions 
for thresholds conditions can be chosen arbitrarily with values 
in the domain of the system. 

Remark 3.1. Our preliminary analysis shows that the dynam- 
ics in a neighborhood of v* can be very complex. We do not 
exclude that a chaos appearance can be observed, and trajecto- 
ries may belong to a fractal, if parameters are not small. It does 
not contradict the zero Lebesgue measure of non-synchronized 
points. Possibly, analysis of non-identical oscillators with not 
small parameters is of significant interests to explore arrhyth- 
mias, earthquakes, chaotic flashing of fireflies, etc. 

Remark 3.2. The time of synchronization for a given initial 
point does not increase if number of oscillators increases (but 
the parameters needed to be closer to zero). This property, pos- 
sibly, can be accepted as a small-world phenomenon. 

4. The simulation result 

To demonstrate our main result numerically, let us consider 
a model of 100 oscillators, which initial values are randomly 
uniform distributed in [0, 1]. Their differential equations are of 
form 

x'j = (4.1+0.01*sort(rand(l,n))-(3.2+0.01*sort(rand(l,n))xj, 

and thresholds 

1 +0.005* sort{rand(\, n)),i= 1,2,..., 100, 

where deviations of coefficients the threshold are also uni- 
formly random in [0,1]. We place the result of simulation with 
e = 0.06 and t = 0.002 in Figure [3] where the state of the sys- 
tem is shown at the initial moment, before the 183-th jump, 
before the 366-th jump and the last is before the 549-th jump. 
That is, it is obvious that eventually all oscillators fire in unison. 



5. Conclusion 

The cardiac pacemaker model of identical and non-identical 
oscillators with delayed pulse-couplings is investigated in the 
paper. We apply the method of investigation proposed in [15], 
which is based on a specially defined map. The map is, in fact 
a Poincare map if one considers the identity of oscillators. Suf- 
ficient conditions are found such that delay involvement in the 
Peskin's model does not change the synchronization result for 
identical and non-identical oscillators JHH, 15\. The result has 
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50 

oscillators 



Figure 3: The state of the model before the first, the 183-th, 366-th, and 
549-th jump is seen. The flat fragments of the graph are groups of oscillators 
firing in unison. 



a biological sense, since retardation is often presents in biolog- 
ical processes and if one proves that a phenomenon preserves 
even with delays, that makes us more confident that the model 
is adequate to the reality. Moreover, the method of treatment of 
models with delay can be useful for neural networks and earth- 
quake faults Jj, S |29l Eo, EjJ analysis. All proved assertions 
can be easily specified with t = 0, to obtain synchronization 
of the Peskin's model for identical QjJ, |2J] and nonidentical 11511 
oscillators. In particular, from the theorem of Section|2]one can 
obtain the result of Example 2.2 of lfl5ll . 
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